rm(list = ls())

library(data.table)
library(estimatr)
library(ggplot2)

load('./data/panel_month_dummies.RData')

## parallel-trends/event study graphs
pre <- panel[month<'2020-01-01']

effect <- vector(mode = 'list', length = length(unique(pre$month)))
ci_low <- vector(mode = 'list', length = length(unique(pre$month)))
ci_high <- vector(mode = 'list', length = length(unique(pre$month)))
month <- vector(mode = 'list', length = length(unique(pre$month)))

for (i in seq_along(unique(pre$month))) {
  print(unique(pre$month)[i])
  
  pre[,post:=ifelse(month==unique(pre$month)[i],1,0)]
  pre[,treat:=post*excess_deaths_jan_d]
  
  out <- lm_robust(hc_pc_asians ~ treat, data=pre, se_type='stata',
                   clusters = panelvar, fixed_effects = ~ panelvar + month)
  
  effect[i] <- out$coefficients
  ci_low[i] <- out$conf.low
  ci_high[i] <- out$conf.high
  month[[i]] <- unique(pre$month)[i]
}

coeffs_pre <- data.frame(effect=unlist(effect), ci_low=unlist(ci_low), ci_high=unlist(ci_high), month=do.call("c", month))

# coefficients post-covid
mon <- seq(as.Date("2020-02-01"), by='month', length.out=2)

effectp <- vector(mode = 'list', length = length(mon))
ci_lowp <- vector(mode = 'list', length = length(mon))
ci_highp <- vector(mode = 'list', length = length(mon))
monthp <- vector(mode = 'list', length = length(mon))

for (m in seq_along(mon)) {
  print(mon[m])
  ppost <- panel[(month<'2020-01-01' | month==mon[m])]
  ppost[,post:=ifelse(month==mon[m],1,0)]
  ppost[,treat:=post*excess_deaths_jan_d]
  
  out <- lm_robust(hc_pc_asians ~ treat, data=ppost, se_type='stata',
                   clusters = panelvar, fixed_effects = ~ panelvar + month)
  
  effectp[m] <- out$coefficients
  ci_lowp[m] <- out$conf.low
  ci_highp[m] <- out$conf.high
  monthp[[m]] <- mon[m]
  
}

coeffs_post <- data.frame(effect=unlist(effectp), ci_low=unlist(ci_lowp), ci_high=unlist(ci_highp), month=do.call("c", monthp))
coeffs <- rbind(coeffs_pre, coeffs_post)
save(coeffs, file='./output/coeff_deaths_jan_backup.RData')